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' Abstract 

' A primary unresolved issue for cosmofogical singularities is whether 

, or not their behavior is locally of the Mixmaster type (as conjectured 

^-H ' by Belinskii, Khalatnikov, and Lifshitz (BKL)). The Mixmaster dy- 

lO . namics first appears in spatially homogeneous cosmologies of Bianchi 

Types VIII and IX. A multiple of the spatial scalar curvature acts as 
O , a closed potential leading, in the evolution toward the singularity (say 

O^' r — > oo), to an (almost certainly) infinite sequence of bounces whose 

$_H . parameters exhibit the sensitivity to initial conditions usually associ- 

ated with chaos. Other homogeneous cosmologies are characterized by 
^ , open (or no) potentials leading to a last bounce as t ^ oo . Such 

models are called asymptotically velocity term dominated (AVTD). 
Here we shall describe a numerical approach to address the BKL con- 
' jecture. Starting with a symplectic numerical method ideally suited 

to this problem, we shall consider application of the method to three 
models of increasing complexity. The first application is to spatially 
homogeneous (vacuum) Mixmaster cosmologies where we compare the 
symplectic ODE solver to a Runge-Kutta one. The second application 
is to the (plane symmetric, vacuum) Gowdy universe on x R. The 
dynamical degrees of freedom satisfy nonlinearly coupled PDE's in one 
spatial dimension and time. We demonstrate support for conjectured 
AVTD behavior for this model and explain its observed nonlinear small 
scale spatial structure. Finally, we study U{1) symmetric, vacuum cos- 
mologies on T'^ X R. These are the simplest spatially inhomogeneous 
universes in which local Mixmaster dynamics is allowed. The Gowdy 
code is easily generalized to this model, although, the spatial differenc- 
ing needed in the symplectic method is not trivial. For AVTD models, 
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we expect the potential-like term in the Hamiltonian constraint to van- 
ish as r ^ oo while in local Mixmaster it becomes (locally) large from 
time to time. We show how the potential behaves for a variety of 
generic U{1) models. 

1 Introduction 

In these lectures, I propose to discuss the application of symplectic numer- 
ical methods [|l|, |2| to the investigation of cosmological singularities. For a 
system whose evolution can be described by a Hamiltonian, the symplectic 
approach splits the Hamiltonian into kinetic and potential subhamiltonians. 
If the subhamiltonians are exactly solvable, these solutions can be used to 
evolve the system from one time to the next. Fortunately, an appropriate 
choice of variables in the standard 3-1-1 Hamiltonian formulation of general 
relativity enables Einstein's equations to be derived from a Hamiltonian ap- 
propriate for the symplectic algorithm (SA). So far, the primary application 
of SA to general relativity has been to determine the asymptotic singularity 
behavior of cosmological models (but see also |^). The SA is well-suited to 
this problem because it becomes exact if the asymptotic behavior is asymp- 
totically velocity term dominated (AVTD) — i.e. the kinetic subhamiltonian 
asymptotically determines the dynamics. 

To study the application of SA to general relativity, we shall consider a 
sequence of models whose variables depend on 0, 1, and 2 spatial dimensions. 
The first case we shall consider is the spatially homogeneous Mixmaster 
universe. (For convenience, we shall consider only the diagonal Bianchi IX 
vacuum model.) Einstein's equations can be obtained by variation of the 
Hamitonian 

2n = -pl + pl + pl + u{n,(5+,(3.) (1) 

where 

-2e'^^+ - 2e-2(/3++v^'3-) - 2e-2(/3+-v^/5-)) . (2) 

Here rt{t) is the logarithm of the cosmological scale factor and measures 
isotropic expansion, while f3± measure orthogonal anisotropic shears with 
and p± respectively canonically conjugate to and (3± The potential 
U is proportional to the spatial scalar curvature {U = e^^ ^i?) and is shown 
in Fig. 1^. Eq. is itself the Hamiltonian constraint W = 0. The properties 
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Figure 1: A typical Mixmaster trajectory in the anisotropy plane with 
horizontal axis /?+ and vertical axis /?_ which are centered on zero. Time 
increases (and the singularity is approached) to the right and downward. 
The number on each frame indicates the axes' scales. The letters label the 
Kasner epochs. The Mixmaster minisuperspace equipotentials are shown 
shaded in gray. 



of this model have been known (more or less) since the late sixties ^. 
The first three terms in U describe triangular walls. As the model evolves 
toward the singularity at = oo, the walls become exponentially steep. 
Within the walls, the system behaves almost as a free particle {U ~ 0). In 
the approach to the singularity, Q itself may be used as the time variable. 
As — > — oo, a fixed value of U, say Uq, moves outward in the P±-plane at 
a speed 1/2 that of the system point 0]. Thus the system evolves with an 
infinite sequence of bounces off the potential. A typical trajectory is shown 
in Fig. 1. While there is no exact solution, each straight line segment can 
be parametrized and a map (the BKL map) derived to link one segment to 
the next i, |, §. There is a long history of numerical simulations ||lO|, [0] 
trying to assess the validity of the BKL map as a descriptor of the dynamics. 
Part of this interest can be traced to the fact that a bounce which leaves 
one of the 120° corners of the potential to move to another corner exhibits 
the sensitivity to initial conditions usually associated with chaos. Whether 
or not Mixmaster dynamics is chaotic remains a topic for both analytic 
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and numerical study [^]. We note that with U = 0, the solution is that 
of Kasner |12|. A solution which is asymptotically Kasner is AVTD. The 
Mixmaster solution is the antithesis of AVTD since there is (presumably) 
no last bounce. The effect of the spatial derivatives that generate U, though 
almost absent during each Kasner epoch, always recurs. 

Our remaining examples are spatially inhomogeneous cosmologies (still 
vacuum for convenience). Long ago, BKL conjectured that the singularity 
in spatially inhomogeneous cosmologies is locally of the Mixmaster type [^. 
Analytic verification of the BKL conjecture has bogged down on the issue 
of setting up the local Mixmaster behavior in a global way. For this reason, 
a numerical approach may be useful. 

While our second model's plane symmetry precludes local asymptotic 
Mixmaster dynamics, it does serve as an excellent laboratory for the SA. 



The Gowdy model on x i? is described by the metric [13| 

+e-^ [e^da^ + 2e^Q da d6 + {e^Q^ + e~^) d6^] (3) 

where A, P, Q are functions of 9, r. We impose spatial topology by 
requiring < 6, a, 6 < 27r and the metric functions to be periodic in 9. 
If we assume P and Q to be small, we find them to be respectively the 
amplitudes of the + and x polarizations of the gravitational waves with 
A describing the background in which they propagate. The time variable r 
measures the area in the symmetry plane with r = oo a curvature singularity. 
Einstein's equations split into two groups. The first is nonlinear ly coupled 
wave equations for P and Q (where ,a = d/da): 

P,rr-e'^^P,ee-e^'' {Q,l-e-^^Q,j) = 0, (4) 
Q,rr-e-^^Q,ee + 2 (P,rQ,r-e-^^P,eQ,e) = 0. (5) 



The second contains the Hamiltonian and ^-momentum constraints which 
can be expressed as first order equations for A in terms of P and Q: 



X,r - [P,l + e-2-p,2 + e2^(Q,2 + Q J )] = 0, (6) 



\,e-2iP,eP,r + e''^Q,eQ,r) = 0. (7) 

This break into dynamical and constraint equations removes two of the 
most problematical areas of numerical relativity from this model: (1) The 
normally difficult initial value problem becomes trivial since P, Q and their 
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first time derivatives may be specified arbitrarily (as long as the total 9 
momentum in the waves vanishes). (2) The constraints, while guaranteed 
to be preserved in an analytic evolution by the Bianchi identities, are not 
automatically preserved in a numerical evolution with Einstein's equations 
in differenced form. However, in the Gowdy model, the constraints are 
trivial since A may be constructed from the numerically determined P and 
Q. For the special case of the polarized Gowdy model (Q = 0), P satisfies a 
linear wave equation whose exact solution is well-known For this case, 
it has been proven that the singularity is AVTD [15|. This has also been 
conjectured to be true for generic Gowdy models ||T^. We shall show in 
Section 4 how the SA applied to this model provides strong support for this 
conjecture. 

Our final model generalizes the plane symmetry to a U(l) symmetry 
[17|. The details of this model will be given in Section 5. U{\) models allow 
local Mixmaster dynamics and thus can be used to test the BKL conjecture. 
In fact, it is possible that any type of allowed cosmological singularity will 
already appear in the U (1) models. The extension of the Gowdy SA methods 
to the U{1) case is straightforward. 



2 Symplectic Methods 

Consider the time evolution of a set of variables X from ti to t2- We can de- 
fine an evolution operator U(t2,ti) such that if f t2 — = is infinitesimal, 
then U must have the form 

U{At)X = + At^^ X. (8) 

But dX/dt = {ff, X} where {i^, X} is the Poisson bracket with H the 
Hamiltonian. Thus U{At) = 1 -|- At{H, } with the empty slot in the 
operator to act on X. In the standard way (by dividing At into n intervals 
and applying Z//(At/n) n times) we obtain the exponentiated form (for finite 
At) 

iY(At) = e^*^^' > = e^*^. (9) 
Suppose H = H1 + H2. Then U = e^*(^i+^2). Consider 

^2)(At) = e^AAt/2)^A,At^MAt/2)_ ^^q) 

Straightforward multiplication shows that the right hand side of (|^) is a 
second order accurate approximation to the evolution operator e^*"^. One 
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evolves X from t to t + At by first evolving with e*-^*/^-*"^^ from t to t + ^ At, 
taking that result and evolving with e^*^^ from t to t + At, and, finally, 
evolving that result with e^^'/^^^^ from t + ^At to t + At. 

Suzuki has shown how to obtain a representation of lA to arbitrary order 
[18|. For example, 

(At) = ^2) (sAt)^2) [(1 - 2s) At]^2) (sAt) (11) 
where s = (2 — 2^/^)^^. In general, 

U{2m) (At) = ^2m-2) (^m At)Z^(2„_2) [(1 " 2Sm) At]^/(2„_2) (^mAt) (12) 

where = (2 — 2-'^/(^™'~-^))~-'^. As a concrete example, consider the Hamil- 
tonian 

H = Hi + H2 = \p^ + V{q) (13) 

where V is an arbitrary potential. Note that both Hi and H2 yield exact 
solutions. For Hi^ 

q{t + At) = q{t)+p{t)At, 

p{t + At) = p{t), (14) 



while for H2, 



qit + At) = qit), 

rlV 

At. (15) 



p{t + At) = p{t) - — 



Since q is constant for the Hamiltonian H2, the solution for p is exact no 
matter how complicated the potential V. Thus the exact solutions ( p^ ) 
and (|l^) are used to evolve from t to t + At according to the prescription 
^/(2)(At). To go to higher order, nothing new is required. The time intervals 
are selected according to the prescription (12), but the same exact solution 
is used. 

Extension to fields q{x,t), p{x,t) is straightforward. With computation 
in mind, define qj = q{xi,P), etc. where Xj = iAx and P = jAt. Then the 
exact solutions are for Hi 

{qi''\pi''') = {ql+piAt,pi) (16) 

and for H2 



6V 



J _ 

5q 




(17) 
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where 6V/6q is the appropriate functional derivative. Again, the exact solu- 
tion exists for H2 no matter how complicated the potential. Of course, one 
must represent the spatial derivative 5V/5q as accurately as possible. 

To represent spatial derivatives to the desired order, evaluate (e.g. in ID) 
the Taylor series for an arbitrary function /(x). Say the 4th order accurate 
expression for d? f /dx^ is desired. Demand that 

ai[f{x + e) + f{x-e)] + 03 [/(x + 2e) + /(x - 2e)] 

= + (18) 

We find oi = |, a2 = — A similar procedure can be used for any term 
in the Taylor expansion. It can also be extended to two spatial dimensions 
where the same coefficients are found. Extension to higher order requires 
f{x + ne) for some n> 2 (depending on order). 

The only remaining issue in the evaluation of {dV /dq)l is whether to vary 
V analytically and then difference the variation or difference V and then vary 
the differenced form. We have chosen the latter. (For complicated V , the 
two are equivalent only to some order.) As a simple example, again consider 

to arise from the variation of = ^ (^) • Since the Hamiltonian 
for this case is H2 = J V\q{x,ty\ dx, the differenced form is (for 2nd order 
accuracy) 

^-ig^W^- '''' 

(The negative of) the variation with respect to (which appears in two 
terms of the sum) yields the standard differenced second derivative (/j+i + 
- 2/,)/(Ax)2. 



3 Mixmaster Model 

Although I have done several computational projects with the Mixmaster 
model, (see references in |^ 20|), none of these have used the SA. Here 



we shall use Mixmaster primarily as an application of SA rather than to 



evaluate Mixmaster parameters or to compute Lyapunov exponents |19|, P0f| . 
The Hamiltonian (||) clearly is in the correct form with Hi = —p'^+p'^ +p'^ 
and H2 = U{P±,Q,) having obvious exact solutions. 

(The variables chosen are not the only ones that have been used. For ex- 
ample, the ADM reduction can be performed by choosing to be the time 
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variable and solving the Hamiltonian constraint (|T|), = 0, for pQ. The 
remaining degrees of freedom for P± are dynamical. Since the constraint is 
automatically preserved, one might argue that this choice is better. How- 
ever, the ability to monitor H as an indicator of the accuracy and stability 
of the solution outweighs the disadvantage of the extra degree of freedom.) 

There is actually some debate over the optimal choice of ODE solver (see 
e.g. [|2^). The choice will involve trade-offs between accuracy and computer 
time. Accuracy is especially important for models with sensitivity to initial 
conditions since a small numerical error could qualitatively change the solu- 
tion. (Sensitivity to initial conditions means that an arbitrarily small change 
in initial data can generate qualitatively different trajectories. Numerical er- 
ror could simulate change in initial conditions.) Here we shall consider only 
comparisons between a 4th order Runge-Kutta algorithm (RKA) and 
4th and 6th order SA's. In all cases, the same initial data set (freely specify 
Q, P±, and p± and solve 7i = ioi pq) is run for up to 2000 time steps. 
An adaptive step size algorithm has been borrowed from [22| for use in all 
these codes. A modification to limit the maximum step size is required to 
avoid seriously over-running a Mixmaster wall. When this eventually oc- 
curred anyway, the computation was stopped. Fig. ^ displays the results 
for the three codes at a single Mixmaster bounce. Note that the 4th order 
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Figure 2: Close up of a Mixmaster bounce in the anisotropy plane. The 
same input file was evolved with a 4th order RKA (open circles) and two SA's 
of 4th (open triangles) and 6th (filled squares) order. The same adaptive 
step size subroutine was used in all cases. 
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SA and RKA have similar performance while the 6th order SA equally well 
represents the solution with many fewer time steps. The 6th order SA easily 
reached r > 10^, more than two orders of magnitude greater than RKA in 
essentially the same number of steps. The much larger time step possible 
with the 6th order SA easily overcomes tje extra computational time needed 
to take a 6th order step as three 4th order ones. It is only now that SA's 
are being applied to the study of Mixmaster |^3| and related models [^. 



4 Gowdy Model on x R 



The SA can be applied to the Gowdy cosmology because the wave equations 
(^) can be obtained by variation of the Hamiltonian 



H 



2tt 





27r 



+ ^ I de 





2 I 2P^ 2 

+ e Q,g 



Hi + H2. 



(20) 



Again equations obtained from the separate variations of Hi and H2 are 
exactly solvable. Variation of Hi yields terms in (Q) containing time deriva- 
tives. These have the exact (AVTD) solution 



P 

Q 

TTp 



-pT + ln[a (1 + C^e^^^)] /3t as t ^ 00, 



a (1 + Ce^f^' 



(1 _ ^2e2/3T 

-2a/3C 



— + C ^Qoasr^cx), 
-> /3 as r — > 00 , 



(21) 
(22) 

(23) 
(24) 



in terms of four constants a, /?, C, and To employ the AVTD solution in 
the SA, the values of P, Q, vrp, and ttq at are used to find q, (3, C, and 
^. These are substuted in ( p[ ) to evolve to new values at according to 
(|To|). Note that evolution with Hi is purely local since there are no spatial 
derivatives. This is advantageous for parallel processing. 

Evolution with H2 is still easy because P and Q are constant. The neces- 
sary differencing has already been discussed in Section 2. For completeness, 
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we give the (2nd order) differenced form of H2 as 



i=0 



The presence of only points and nearest neighbors in (25) also yields easy 
parallelization. The exponential prefactor e^^"^ in H2 makes plausible the 
conjectured AVTD singularity. However, P ^ vt (for u > 0) as r — > 00, 
(where from (^) v = (5). If v > 1, the term e~'^'^ e^^ Q in ( pOD can grow 
rather than decay as r — > 00. This has led to the conjecture that the AVTD 
limit requires v < 1 everywhere except, perhaps, at a set of measure zero 
(isolated values of 0) 111 



Implementation of the SA is completely straightforward. The published 



results [24] are based on a code which is 4th order accurate in both time 
and space. Actually, accurate results rely most strongly on accurate repre- 
sentation of the spatial derivatives. The time step must be chosen to satisfy 
a Courant condition. To test the code, we note that it is possible to trans- 
form an exact solution Pq = ^0 (e~^) cosO of the polarized case (for 1^(2;) 
an irregular Bessel function of order n) into the pseudo-unpolarized solution 
P = In cosh Pq) Q = tanhPo which satisfies the full equations (^). Substi- 
tution shows that the nonlinear terms (which must be absent in a polarized 
model) miraculously cancel out. However, the code is unaware of this a 
priori and all terms are present. The results are published elsewhere [^] 
and demonstrate the need for the 4th order code. 



Most of the remaining results [24, ^] use the initial data P = 0, ttp 



vqcosO, Q = cos 9, and ttq = 0. This model is actually generic for the 
following reasons: The cos 9 dependence is the smoothest nontrivial possi- 
bility. With cos n9, the solution is repeated n times on the grid yielding the 
same result with poorer resolution. The amplitude of Q is irrelevant since 
the Hamiltonian (|20|) is invariant under Q — > pQ, P ^ P — \n p. This also 
means that any unpolarized model is qualitatively different from a polarized 
{Q = 0) one no matter how small Q. 

The accuracy and stability of the code easilty allow verification of the 



conjectured AVTD behavior [24|. A plot of the maximum value of v vs r 
(Fig. ^) shows strong support for the conjecture that v < 1 in the AVTD 
regime. However, Fig. ^ also shows that a simulation at higher spatial 
resolution begins to diverge. Normally, failure of convergence signals nu- 
merical problems. Here something different is occurring. The evolution 
of spatial structure in P depends on competition between the two nonlin- 
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Figure 3: Plot of Vmax vs r. The maximum value of v is found for two 
simulations with 3200 (solid line) and 20000 spatial grid points (circles) 
respectively. The horizontal line indicates v = 1. 



ear terms in the P wave equation. Approximating the wave equation by 
P^^^ +either of the nonlinear terms = 0, a first integral can be obtained. 
The two potentials are Vi = Tr^e"^^ and V2 = e~^'^e^^(5j- Non-generic 
behavior can arise at isolated points where either Q,g or ttq vanishes. Say 
such a point is ^o- The finer the grid, the closer will be some grid point 
to ^0- Thus non-generic behavior will become more visible on a finer grid. 
Detailed examination shows that the differences seen in Fig. |3| are due to 
the slower decay of u to a value below unity at isolated grid points in the 
higher resolution simulation. 

Details of the high resolution simulation (for < 6 < 27r) are shown in 
Fig. ^. The narrow peaks in P occur where Q,e^ 0. Generically, if P ~ u r 
and V > 1, the potential V2 dominates. The relevant first integral of (^) is 

^j' + e2^Qj= const (26) 

where Z = P — t. A bounce off V2 yields dZ/dr — > —dZ/dr 01 v — 1 ^ 1 — v. 
If the new v < 0, then Vi dominates yielding the first integral 

P,l +e-2^7r^ = const (27) 

causing P,^ —P,t 01 v ^ —v. Eventually, bouncing between potentials 
gives V < 1. However, if Q,g 0, but is not precisely zero, it takes a long 
time for the bounce off V2 to occur. Precisely at ^0 (where Q,o = 0), v > 1 
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Figure 4: P (dashed line) and Q (solid line) vs6 at t = 12.4 for the standard 
initial data set with vq = 5 for < 9 < ir for a simulation containing 20000 
spatial grid points in the interval [0,27r]. The numbers on the graph refer 
to the most interesting features. Peaks 1, 2, and 3 in P are essentially the 
same in that they occur where Q,g ~ 0. of Q. Peak 4 shows an apparent 
discontinuity in Q where ttq « 0. 




persists. The apparent discontinuity in Q is not a numerical artifact. It 
occurs where P < and ttq ^ 0. Since Q,r = e~^^7rQ and ttq « c{6 — 6q) (if 
ttq = at 6*0), Q,T grows exponentially in opposite directions about ^o- The 
potential Vi drives P to positive values unless ttq = 0. Thus this feature 
will narrow as the simulation proceeds. 

Finally, we report a strange, and yet not understood, scaling of spatial 
structure in P with the parameter vq in the initial data. In our standard 
initial data set, greater values of vq lead to the appearance of additional 
spatial structure in a shorter time. The rate of structure formation decreases 
and then stops as AVTD is approached. One may count the number of 
peaks in P (a peak is crossed if {Pi+i — Pi){Pi — Pi-i) < and Pi+i < Pi) 
during the simulation. The scaling is best for the time at which the 5th 
peak appears although it is also seen for T3 and ry (the even nature of the 
solution causes two peaks to appear at once except at = vr). A plot of 
l/rs vs vq yields a straight line as shown in Fig. ^ which may be described 
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as Ts = a{vo — Vq^) ^ where, if vq = Vq', the 5th peak does not appear until 
r = oo. Even more surprisingly, almost the same line is obtained for initial 




Figure 5: Scaling in the Gowdy model. Plot of l/rs, the inverse of the 
time at which the 5th peak appears in P vs vq. Two cases are shown for 
initial data P = 0, Q = cos9: (1) vrp = {vq / \/2) cos 6 , ttq = {vq/^/2)cos6 
is indicated by filled circles; (2) vrp = vqcosO, ttq = is indicated by open 
squares. 

data Tip = {vq/\/2) cos 9, ttq = (vq/V^) cos 9 rather than our standard data. 
Explanation of this scaling is still in progress. 



5 U{1) models 



Moncrief has shown |17] that cosmological models onT^ x R with a spatial 
U{1) symmetry can be described by five degrees of freedom {x,z,A,{p,uj} 
and their respective conjugate momenta {px,Pz,PA,P,f}- All variables are 
functions of spatial variables u, v and time r (related to distance in the 
symmetry direction). If we define a conformal metric gab in the u-v plane 
as Qab = e^eabix, z) where 



/ + e-2^(i + a:)2 e^' + e-^'ix^ - 1) 
^ g2. ^ e-2^(x2 - 1) + e-2-(l - xf 



(28) 
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has unit determinant and choose the lapse N = e^, Einstein's equations can 
be obtained by variation of 

= H1+H2. (29) 

Note particularly that 

Hi = H^{-2z,x) + H^i-2ip,u;) + H[{A) (30) 

where H^{P,Q) is the kinetic part of the Gowdy Hamiltonian ([20|). Of 
course, Hf is just a free particle Hamiltonian for the degree of freedom 
associated with A. This means that not only are the equations from Hi 
exactly solvable but also that the Gowdy coding can be used with essentially 
no change. The potential term H2 is very complicated. However, it still 
contains no momenta so its equations are trivially exactly solvable. Thus, 
at least in principle, the extension of the Gowdy code to the two spatial 
dimensions of the U{1) code is completely straightforward. 

There are three complications, however, which cause the U{1) problem 
to be more difficult. The first involves the initial value problem (IVP) — the 
constraints must be satisfied on the initial spacelike slice. The constraints 
are 

no = n-2pA = (31) 
(where 7i is the density in (p9|)) and 

Ha = -2TXa± + PkKa -PA,a +PV,a +ru;,a = (32) 

where vf^ is in the 2-space with metric Cab and is linear in px and pz with 
each term containing one or the other. Moncrief has proposed a particular 
solution to the IVP. First, identically satisfy TC = hy choosing 

Px = Pz = V>,a = ^^,a= ; pA = ce^ (33) 

for c a constant. Then, solve Tl^ = for either r or p. Solution is possible 
for c > Cmin such that or p"^ > 0. This allows x, z, A, and p or r to be 
freely specified. (Without loss of generality, it is possible to set x = z = 
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initially to yield Cab flat. Such a condition may always be imposed at one 
time by rescaling u and v.) 

The second difficulty also involves the constraints. While the Bianchi 
identities guarantee the preservation of the constraints by the Einstein evo- 
lution equations, there is no such guarantee for differenced evolution equa- 
tions. At this stage of the project, we monitor the maximum value of the 
constraints vs r over the spatial grid but do nothing else to try to stay on 
the constraint hypersurface. 

The third difficulty, and the one that is proving to be the greatest ob- 
stacle, is instability associated with spatial differencing in two dimensions. 
The differencing algorithm of section 2 for sixth order accurate expressions 
(based on second order ones) is used. The difficulties we encounter are ap- 
parently common to nonlinear equations in two or more dimensions. In an 
attempt to control the instability, we have introduced a form of spatial av- 
eraging. At the end of every time step, each value of every variable, say , 
is replaced by 



This expression gives ^ij = + 0(A^). In both test cases and generic mod- 
els, the averaging procedure has allowed the code to run longer. However, 
the fact that 7^ ^ij can lead to deviations of the numerical solution from 
the correct one. Fortunately, by comparing runs with and without averag- 
ing, these artifacts are easy to identify. We shall see some examples of how 
averaging can allow the code to run long enough for a conclusion about the 
asymptotic singularity behavior to be drawn. 

Moncrief has provided a test case for the U (1) code. It again starts with 
a polarized Gowdy solution and transforms it as either a ID — > w or v) 
or 2D {9 f{u,v)) test problem to satisfy the U{1) equations (including 
the constraints). As a ID example, the agreement is excellent and the code 
can be run to large r. Difficulties arise in the 2D test problem in regions 
where the spatial derivatives are large. In application of the U (1) code to 
generic models, AVTD models can display nonlinear wave interactions before 
settling down to U 0, z,ip,A ^ const r, and x,a; — const. Increasing 
spatial resolution will yield narrower (and steeper) structures and thus may 
not help to cure instabilities due to steep gradients. 

Nonetheless, conclusions can be drawn for generic models in our re- 
stricted class of initial data. We shall consider the models as representative 



-| (Ci+lj+1 + + + • 



(34) 
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of subclasses of the data. Models with r = uj = are called polarized. This 
condition is compatible with the above solution and is preserved (identically) 
by the (analytic and numerical) evolution equations. Grubisic and Moncrief 
have conjectured that these polarized models are AVTD |26]. Therefore, 
the first model is chosen to be polarized. It exhibits the conjectured AVTD 
behavior as shown in Figs. |6| and 0. The second is generic with p given 





Figure 6: Frames of U{u,v,t) for the polarized model x = z = A = 
sin « sin?;, = 12e^, u> = r = 0. Time increases to the right and down- 
ward. The final frame corresponds to [/ « everywhere. 

and the Hamiltonian constraint solved for r in the IVP. Averaging as in 
( p^ allows the model to be followed to the point where only artifacts have 
U ^ 0. This is shown in Fig. |8|. The last model has r given with p obtained 
by solving the Hamiltonian constraint in the IVP. Models of this type are 
less stable, probably due to the growth of a steep feature in to which does 
not appear in the other cases. For this reason, the parameters must be kept 
small. Fig. ^ shows that C/ — > except near artifacts where no statement 
can be made. 
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Half-Polarized Model 




Figure 7: Surface plot of all U (1) variables for the line u = —v vs r for 
the polarized model. Note agreement with predicted AVTD behavior that 
X becomes constant in r while A, 99, and z grow linearly. 

6 Conclusions 

The application of the SA to Einstein's equations for collapsing universes 
allows the nature of the singularity to be studied. While no real attack 
on Mixmaster with SA has been made, the method offers the potential for 
efficient, accurate evolution of this model. Application of SA to the Gowdy 
model has yielded strong support for its conjectured AVTD singularity and 
has allowed the discovery and study of interesting small scale spatial struc- 
ture and scaling. 

Further progress in understanding the generic singularity of U{1) cos- 
mologies requires imporvements in handling steep spatial gradients. Never- 
theless there is strong support that (at least within our restricted class of 
initial data) polarized models are AVTD. There is also support for AVTD 
behavior in all generic models studied so far. Mixmaster-like bounces have 
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Figure 8: Frames of U{u,v,t) for the generic model x = z = 
cos cost;, A = sin u sin t;, pA = 14e^, p = 10 cos « cost; with averaging. 



not been seen. (The activity due to nonlinear wave interaction seen early in 
the simulations is similar to that in Gowdy models.) Several factors could 
account for this: (1) the BKL conjecture is false; (2) the simulations have 
not run long enough; (3) Mixmaster behavior is present but hidden in our 
variables; or (4) our class of initial data is insufficiently generic. All these 
possibilities will be explored in studies in progress. 
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Figure 9: Frames of U{u, v, r) for generic model x = z = 0, A = 
.Icoswcosf, PA = 2.1e^, r = cos « cos The diagonal features in the 
final frames are numerical artifacts. 



Appendix 

This Appendix contains a FORTRAN subroutine which uses the 6th order 
SA to 

integrate Einstein's equations for the Mixmaster universe: 

dn df3± 



where 



_2e^^+ - 2e-2(/3++^/5-) - 2e-2(/3++V3/3-). (37) 
subrout ine sa6 (x , t , tau , param , nstate , xout ) 
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*TAKES ONE 6TH ORDER SA STEP FROM T TO T + TAU 

* 

*VARIABLES: X(l) = BP, X(2) = BM, X(3) = W 
*X(4) = PP, X(5) = PM, X(6) = PW 

* 

*XOUT CONTAINS THE OUTPUT VALUES OF THE VARIABLES 
*PARAM = SQRT(3) , NSTATE = 6 

* 

*XTEMP AND XS ARE DUMMY ARRAYS 

*V IS THE POTENTIAL, DVDBP,DVDBM ARE THE GRADIENTS 
*TERM_I IS AN EXPONENTIAL 

* 

*S, SI ARE THE 4TH,6TH ORDER SUZUKI PARAMETERS 
*DELT,DELT1 STORE T INTERVALS 
parameter (mstate=20) 

real*8 x(6) ,param(l) ,xout (6) ,xtemp(6) ,xs(6) 
real*8 terml , t erm2 , termS , term4 , termS , term6 , 
&v , dvdbp , dvdbm , a , s , delt (3) , si , delt 1 (3) 

* 

*DOUBLE PRECISION IS NECESSARY 
*S WAS COMPUTED WITH MATHEMATICA 

3=1.35120719195965 

sl=l . dO/ (2 . dO-2 . dO** ( . 2d0) ) 

* 

*STORE THE VARIABLES IN A DUMMY ARRAY 
do i=l,6 
xs(i) = x(i) 
enddo 

* 

*ONE 6TH ORDER STEP IS THE PRODUCT OF 3 
*4TH ORDER STEPS 

* 

♦COMPUTE THE INTERVALS FOR THE 4TH ORDER STEPS 

deltl(l) = sl*tau 

deltl(2) = (I.d0-2.d0*sl)*tau 

deltl(3) = sl*tau 
*THE DELTl SUM TO TAU 
* 

*EACH 4TH ORDER STEP IS THE PRODUCT OF 3 
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*2ND ORDER STEPS 
* 

do idell=l,3 

* 

*EACH 4TH ORDER INTERVAL IS DIVIDED INTO THREE 
* INTERVALS FOR 2ND ORDER 

* 

delt(l) = s*deltl(idell) 

delt (2) = (1 . dO-2 . dO*s) *delt 1 (idell) 

delt(3) = s*deltl(idell) 

* 

do idel =1,3 

^PERFORM A SECOND ORDER SA STEP FOR EACH INTERVAL 

* 

half _tau=0 . 5*delt (idel) 
a = param(l) 

* 

*EVOVLE WITH H_l FOR HALF A TIME STEP 

* 

xtemp(l) = xs(l)+xs(4)*half _tau 
xtemp(2) = xs(2)+xs(5)*half _tau 
xtemp(3) = xs(3)-xs(6)*half _tau 

* 

*EVALUATE THE EXPONENTIALS IN THE MOST ACCURATE WAY 
terml=exp (4*xteiiip (3) -8*xtemp ( 1 ) ) 
term2=exp (4*xtemp (3) +4*xtemp ( 1) ) 
terni3=exp (4*xt emp (3) +4*xtenip ( 1 ) +4*a*xtemp (2) ) 
term4=exp (4*xtemp (3) -2*xtenip ( 1 ) -2*a*xtemp (2) ) 
term5=exp (4*xteiiip (3) +4*xtemp ( 1 ) -4*a*xteinp (2) ) 
term6=exp(4*xtemp(3)-2*xtemp(l)+2*a*xtemp(2) ) 

* 

♦COMPUTE V AND ITS GRADIENTS 

v=terinl+terin3+term5-2 . * (term2+term4+term6) 
dvdbp=4 . * (-2 . * (terml+term2) +term3+terni4 
&+term5+term6) 

dvdbm=4 . *a* (term3+term4-term5-term6) 

* 

*EVOLVE WITH H_2 FOR A FULL TIME STEP 
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xs(4) = xs(4) - . 5*dvdbp*delt (idel) 
xs(5) = xs(5) - .5*dvdbm*delt(idel) 
xs(6) = xs(6) - 2.*v*delt(idel) 

* 

*EVOLVE WITH H_l FOR HALF A TIME STEP 
* 

xs(l) = xtempd) + xs(4)*half _tau 
xs(2) = xtemp(2) + xs(5)*half_tau 
xs(3) = xtempO) - xs(6)*half _tau 

* 

enddo 
enddo 

* 

*RECORD THE DUMMY ARRAY 
do i=l,6 
xout(i) = xs(i) 
enddo 

* 

return 
end 

OTHER CODE FEATURES: 

Input w, bp, bm, pp, pm and solve H for pw. 

A positive square root yields a collapsing universe. 

Use an adaptive step size but require tau < .01 t. 
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